OHI British Columbia | OHI Science | Citation policy
OHI British Columbia | OHI Science | Citation policy
knitr::opts_chunk$set(fig.width = 6, fig.height = 4, fig.path = 'Figs/',
echo = TRUE, message = FALSE, warning = FALSE)
library(sf)
dir_git <- '~/github/ohibc'
source(file.path(dir_git, 'src/R/common.R')) ### an OHIBC specific version of common.R
dir_rgn <- file.path(dir_git, 'prep/regions') ### github: general buffer region shapefiles
dir_anx <- file.path(dir_M, 'git-annex/bcprep')
### goal specific folders and info
goal <- 'fis'
scenario <- 'v2017'
dir_goal <- file.path(dir_git, 'prep', goal, scenario)
dir_goal_anx <- file.path(dir_anx, goal, scenario)
dir_spatial <- file.path(dir_git, 'prep/_spatial')
dir_dfo_data <- file.path(dir_anx, '_raw_data/dfo_khunter')
### provenance tracking
library(provRmd); prov_setup()
### support scripts
source(file.path(dir_git, 'src/R/rast_tools.R'))
### raster plotting and analyzing scriptsThis script collects spatial boundaries for RAM stocks, using stock IDs from the RAM prep script and the boundaries created by Chris Free: https://marine.rutgers.edu/~cfree/ram-legacy-stock-boundary-database/
From this, we generate for each stock a region-by-region proportional area to apportion RAM catch to OHIBC regions.
Reference: https://marine.rutgers.edu/~cfree/ram-legacy-stock-boundary-database/
Downloaded: 7/5/2017
Description: Boundaries for RAM stocks
Format: ESRI shapefile format.
For each RAM stock identified in the script 1_data_prep_fis_ram.Rmd, identify the appropriate shapefile for boundaries, either from the Chris Free datasets or from DFO (if there is a difference).
The .zipped set of shapefiles seems to be named according to the RAM assess_id field. Copy the relevant set to git-annex, in the fis/v2017/stock_boundaries/ram folder.
stock_list <- read_csv(file.path(dir_goal, 'output/ram_catch.csv')) %>%
select(stock_id, stock_name, ram_area_id, ram_area_name) %>%
distinct()
ram_spatial_dir <- file.path(dir_anx, '_raw_data/ram_fisheries/d2017/spatial')
boundary_list <- readxl::read_excel(file.path(ram_spatial_dir, 'ramldb_v3.8_stock_boundary_table_v2_formatted.xlsx')) %>%
filter(stockid %in% stock_list$stock_id)
shp_list <- list.files(file.path(ram_spatial_dir, 'ramldb_boundaries'), full.names = TRUE)
boundary_list <- boundary_list %>%
select(assessid, stockid, stocklong, zone_col, zones, notes) %>%
mutate(dir = file.path(ram_spatial_dir, 'ramldb_boundaries'),
shp = assessid,
file_exists = file.exists(file.path(dir, shp)))
bc_bounds_dir <- file.path(dir_goal_anx, 'stock_boundaries/ram')
unlink(bc_bounds_dir, recursive = TRUE)
dir.create(bc_bounds_dir)
y <- lapply(boundary_list$shp, FUN = function(x) {
### x <- boundary_list$shp[1]
stockname <- x %>%
str_replace('-PAC-', '-') %>%
str_split('-') %>%
unlist() %>% .[2]
file.copy(from = file.path(ram_spatial_dir, 'ramldb_boundaries', paste0(x, c('.shp', '.prj', '.shx', '.dbf'))),
to = file.path(bc_bounds_dir, paste0(stockname, c('.shp', '.prj', '.shx', '.dbf'))))
})For each RAM stock, intersect the RAM boundary polygons (from Chris Free) with the OHIBC regions, then calculate area (in km^2) within each OHIBC region and proportional areas, relative to the overall stock polygon area (including stock extents outside the BC EEZ). The output layer for toolbox use will simply contain OHIBC region, stock ID, and stock-in-region area in km^2.
Note that most species do not actually get fished in the Pacific Offshore region - coastal species and/or demersal species. For these, we will remove area associated with the Pacific Offshore region and recalculate total area based on the continental shelf portion of the polygon instead.
Two exceptions: Halibut and Albacore. These will base the regional catch on the total area of the fishery polygon. The catch for Halibut associated with Pacific Offshore will be dropped.
ram_areas_file <- file.path(dir_goal, 'ram', '2_ram_stock_to_ohibc_areas.csv')
reload <- TRUE
if(!file.exists(ram_areas_file) | reload) {
ram_bounds_dir <- file.path(dir_goal_anx, 'stock_boundaries/ram')
shp_list <- list.files(ram_bounds_dir, pattern = '.shp$', full.names = FALSE) %>%
str_replace('.shp$', '')
ohibc_sf <- read_sf(dir_spatial, 'ohibc_rgn') %>%
st_transform(4326)
intsx_ram <- function(ram_layer) {
### ram_layer <- shp_list[2]
ptm <- proc.time()
ram_sf <- read_sf(ram_bounds_dir, ram_layer) %>%
mutate(tot_area_m2 = st_area(geometry))
intsx_sf <- st_intersection(ohibc_sf, ram_sf) %>%
select(rgn_name, rgn_id, assessid, stockid, tot_area_m2)
if(nrow(intsx_sf) > 0) {
intsx_sf <- intsx_sf %>%
mutate(area_m2 = st_area(geometry))
if(!str_detect(ram_layer, 'PHALNPAC|ALBA')) {
### subtract pacific offshore and readjust total area
intsx_sf <- intsx_sf %>%
filter(rgn_id != 7) %>%
mutate(tot_area_m2 = sum(area_m2))
}
if(str_detect(ram_layer, 'PHALNPAC')) {
### subtract pacific offshore, but leave total area untouched
intsx_sf <- intsx_sf %>%
filter(rgn_id != 7)
}
intsx_df <- intsx_sf %>%
as.data.frame() %>%
select(-geometry) %>%
mutate(elapsed = (proc.time() - ptm)[3])
return(intsx_df)
} else {
return((proc.time() - ptm)[3])
}
}
# ram_ohibc_all_list <- vector('list', length = length(shp_list)) %>%
# setNames(shp_list)
# for(ram_layer in shp_list) {
# message(ram_layer)
# ram_ohibc_all_list[[ram_layer]] <- intsx_ram(ram_layer)
# }
ram_ohibc_all_list <- parallel::mclapply(shp_list, FUN = intsx_ram, mc.cores = 18) %>%
setNames(shp_list)
ram_ohibc_dfs <- ram_ohibc_all_list[sapply(ram_ohibc_all_list, class) == 'data.frame']
ram_ohibc_all <- ram_ohibc_dfs %>%
bind_rows() %>%
group_by(stockid) %>%
mutate(area_km2 = area_m2 / 1e6,
a_prop = area_m2 / tot_area_m2,
tot_area_km2 = tot_area_m2 / 1e6) %>%
select(-elapsed, -area_m2, -tot_area_m2)
write_csv(ram_ohibc_all, ram_areas_file)
}ram_ohibc_all <- read_csv(ram_areas_file)
DT::datatable(ram_ohibc_all)ram_ohibc_all <- read_csv(file.path(dir_goal, 'ram', '2_ram_stock_to_ohibc_areas.csv')) %>%
select(rgn_id, stockid, area_km2, a_prop)
write_csv(ram_ohibc_all, file.path(dir_goal, 'output', 'ram_stock_area.csv'))ram_bounds_dir <- file.path(dir_goal_anx, 'stock_boundaries/ram')
shp_list <- list.files(ram_bounds_dir, pattern = '.shp$', full.names = FALSE) %>%
str_replace('.shp$', '')
ohibc_sf <- read_sf(dir_spatial, 'ohibc_rgn') %>%
st_transform(4326)
for(ram_layer in shp_list) {
### ram_layer <- shp_list[10]
message('Generating map for ', ram_layer)
ram_sf <- read_sf(ram_bounds_dir, ram_layer)
ram_map <- ggplot() +
ggtheme_plot() +
geom_sf(data = ohibc_sf, alpha = .5, fill = 'slateblue', color = 'blue', size = .25) +
geom_sf(data = ram_sf, alpha = .3, fill = 'red', color = 'red', size = .25) +
labs(title = ram_layer)
print(ram_map)
}
prov_wrapup(commit_outputs = FALSE)